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Abstract 

The discrete gradient approach is generahzed to yield integral preserving methods for 
differential equations in Lie groups. 

1 Introduction 

Our point of departure is the system of differential equations 

x^F{x)^f{x)-x^x-~f{x), (1) 

where the unknown x = x{t) is a curve on some Lie group x{t) e G and the dot over x signifies 
differentiation with respect to t. The map f ' G ^ q where Q = TeG is the Lie algebra correspond- 
ing to G, and the dot should be interpreted as the derivative of right (resp. left) multiplication, 
e.g. f{x)'X := TeRxf{x). Such equations occur for instance in mechanical systems where the Lie 
group could be the special orthogonal group SO (3) or the special Euclidean group SE{3). We 
shall consider in particular the case where the system ([T]) possesses one or more first integrals. 
Here we define a first integral to be any function H : G which is invariant on solutions 

^^H(x(t)) = {dH,F) = 0. 

Any differential equation ([T]) on a Lie group having H as first integral can be formulated via a 
bivector (dual two-form) uj and the differential of H as 

X = F{x) = uj{dH, ■) = dHjuj (2) 

A Riemannian metric on G defines an inner product {•^•)x on every tangent space T^G which is 
also varying smoothly with x. With such a metric one may define the Riemannian gradient vector 
field as the unique vector field gradi^ satisfying {dH,v)x = {gTdidH\x,v)x for every v e T^G. An 
example of a bivector uj to be used in ([2| is then provided by means of the wedge product 

gradi^ A F 
llgradi^ P 

Note that u is not uniquely defined by F and H. For a choice x = (xi, . . . , x^^) of local coordinates 
on the group, we may write 

d 

^\x= X! ^ij{x)dxi Adxj, di^l^ = X - — dxk (3) 

l<i<j<d k=l 
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In this way, we find the coordinate version of (§ 

x = S^VH (4) 

for the skew- symmetric d x (i-matrix S^j = u;^ - u) and u) is the strictly upper triangular matrix 
with entries ujij{x)^ 1 < i < j < d. The formulation Q has been the starting point for energy 
preserving integrators devised by Gonzalez [8 . In fact, McLachlan et al. [13] showed that under 
relatively general circumstances, any vector field F on with a first integral H can be written 
in the form 

F{x) = S{x)VH (5) 

for some skew-symmetric matrix S(x). 

Note that the bivector is not required to be nondegenerate. For Hamiltonian systems, the 
Hamiltonian itself is a first integral and an accompanying bivector can be inferred from the 
symplectic two-form Q by inversion. In coordinates one can represent the symplectic two-form 
by a skew-symmetric matrix Sq in a similar way as in (|3| and Q with respect to the basis 
{dxi A dxj, i < j}. By definition, this matrix will be invertible, and its inverse is precisely 

The formulation ([2| is easily generalised to the case with k independent first integrals 
i^i, . . . , i^/c. We may now replace the bivector with a A: + 1-vector uj and write ([T]) as 

x = F{x)=uj{dHi,...,dHk, •) 

Also in this case we can find, by means of a Riemannian structure, an example of a feasible 
k + 1-vector 

^ = — 7T77 TTTT where ujq = gradi^i A ••• A gradiJ/e 

cJo(diii,...,dii/e) 

Earlier work on energy-preserving methods on Lie-groups was presented in [12], see also the 
references therein. In this paper we shall generalise the notion of discrete gradient methods to 
the situation where the phase space is a Lie group. The paper extends results presented in [4^ 
by also considering high order methods, a larger class of manifolds and further examples. 



2 Discrete differentials in Lie groups 
2.1 A review of the situation in Euclidean space 

The idea of discrete gradient methods, is to consider some approximation to the exact gradient 
in ([5| xR^ satisfying the following two conditions 

H{v) - H{u) = VH{u, v)^{v - u), Mu, v e (6) 

VH(u,u) = VH(u), ViiGR^ (7) 

Many such discrete gradients have been proposed in the literature, and we give here a few 
examples. The averaged vector field discrete gradient is discussed for instance in [13j and more 
recently in the PDE setting [71 13] 

VH{u,v)= [\H{{l-Ou + ^v)dC (8) 
Jo 

Another discrete gradient which was proposed in [8] is the midpoint gradient 
_ /u + v\ H(v) - H(u) -VH(^f (v-u) 

VH{u, v) = Vh{^)^ ^ ^^^^ ^ {v - u) (9) 
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An integral preserving integrator for ([5| is readily given as 

n+1 _ n 

- — — ^ = 5(x",x"^^)Vi/(x",x"^') (10) 
h 

where is a skew- symmetric matrix which approximates S. It is required to satisfy the consis- 



tency condition S{u,u) = S{u). An easy calculation, using (|6]) and (10) now shows that 

the last identity follows since S is skew-symmetric. The skew-symmetric matrix S used in the 



integrator is not unique, and the freedom can be used to improve the approximation (10) in 
various ways, for instance to increase its order of convergence. 

2.2 The Lie group setting 



For Lie groups, none of the definitions mh or (10) can be used since they both involve vector 



space operations, not generally defined on Lie groups. In a coordinate free setting, we also find 
it more convenient to replace the gradient and its discrete counterpart by dual quantities, the 
differential as indicated in (|2|. Below, we also use the notion of a discrete differential rather 
than a discrete gradient. As is the tradition for Lie group integrators [l0l|5], approximations are 
introduced through some finite dimensional action and a trivialisation principle is applied. By 
this, we mean that tangent vectors at some x g G, i.e. Vx ^T^G can be represented via either left 
or right translation of a vector ^ e T^G = where g is the Lie algebra of G. In the present paper 
we will just for convenience choose right translation. Defining the right multiplication operator 
Rx ' G ^ G^ where Rxy = y - we shall use the notation 

and similarly, any vector Px ^T^G is identified by some /i e through 

{Px^Rx^O = {KPx^O = (m.O i-e. li^RlVx, 

where (•,•} is some duality pairing between T^G and X^G, as well as between and q. For 
example, if the Lie group G and its Lie algebra ^ both are realized as m x m-matrices, the 
dual elements can also be represented as matrices and the duality pairing could be given as 
(p, v) = tTdiCe{p^v). In this case one simply has Vx = Rx^£, = and /i = RxPx = Px'X^ • We also 
note that in Euclidean space (G = and group operation is +), left and right translation are 
both realised by the identity map, e.g. Rx^£, = 

We shall introduce the trivialised discrete differential of a function as a map dH : GxG ^ g"" 
satisfying the following identities generalised from (|6| and ^ 

H{v) - H{u) = {dH{u, v),log{v • u'^)) (11) 
dH(x,x)=RldHx (12) 

By identifying vectors and co- vectors in Euclidean space through the standard inner product, 
and noting that the logarithmic map log : G ^ Q in this setting is simply the identity map, 
thus \og{v • u~^) = V - we recover the standard discrete gradient conditions (|6| and ^ when 
Euclidean space is chosen as our Lie group. We also need to introduce a trivialised approximation 
to the bivector a; in ([2|. For this purpose we define, for any pair of points (u^v) e G x G, an 
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exterior 2-form on the linear space q"" which we denote by uj{u^v) thus u : G x G ^ A'^{q''). We 
impose the consistency condition 

Of course, in practice, Co need only be defined in some suitable neighborhood of the diagonal 
subset X e G}. Introducing coordinates, the form Co plays a similar role as the skew- 

symmetric matrix S{x^ ^x^^^) in (10). We may further define our numerical method as follows 

x^-'^ =exp(/iF(x^,x^-'^))-x^, = di^(x^,x^-'^) JcD(x^,x^-'^) (13) 

For the reader who is unfamiliar with the notation, the definition of F{x^ ^x^'^^) e q using 
coordinates as in Q would be F{x^,x^'^^) = Suj{x'^ , x'^'^^)dH {x'^ , x'^'^^) where Suj{x^ ,x^'^^) is a 
skew-symmetric matrix approximating Suj{x'^) and where we have also expressed di^(x^, x^"^^) 
in coordinates. 

From the defining relation (11), it follows immediately that the method preserves H since 

Hix""^^) - H{x'') = (di7(x",x"^^),log(x"^^(x")-^)) = /i(di7(x",x"^^),F(x^,x^^^)) 
= Cd{x'',x''^^)(dH{x'',x''^^),dH(x'',x''^^)) = 

2.3 Examples of trivialised discrete differentials 

An example of a trivialised discrete differential, generalising (|8| is 

dH{u,v) = £ R;^^^dHe(^^jd^, e{0-exp{^log{vu-'))-u (14) 

Here, we have introduced a straight line in the Lie algebra between the points and a = log{v'U~^). 
By applying exp to each point on the curve and multiplying by li, we obtain a curve on the Lie 
group between the points u and this is Finally the (trivialisation of the) differential dH is 
averaged along this curve to obtain the AVF type of trivialised discrete differential. Considering 
with u and v interchanged yields £{1 -^), and from this it easily follows that dH{u^v) = 
dH{v,u). 

The Gonzalez midpoint gradient can be generalised as well, for instance by introducing an 
inner product on the Lie algebra, we denote it (•,•)• We apply "index lowering" to any element 
7^ e by defining rj^ g to be the unique element satisfying (?7\ C) - ivX) C ^ We can 

then introduce a generalisation of ([9| as 

dH{u,v)=RldH\, + ^^ \^ ; ^ ^^r]\ r] = \og{vu ^), (15) 

where ce G^ is some point typically near u and v. One may for instance choose c = exp (7^/2) • 
which implies symmetry, i.e. dH{u^v) = dH{v^u). 

We have now presented two examples of trivialised discrete differentials which are both sym- 
metric in the two arguments u and v. One observes from the definition ([l3| that the integrator 
itself is symmetric if dH{u^v) = dH{v^u) and uo{u^v) =Cd{v,u) for all pairs (u^v). 



3 More general manifolds 

More interesting examples of mechanical systems can be found for instance in the larger class 
of homogeneous manifolds. There are various ways to devise discrete gradient methods in this 
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setting, or for even more general classes of manifolds. From now on, we assume that M is a 
smooth manifold for which a retraction map is available, retracting the tangent bundle TM into 
M. This is a very basic and straightforward approach, certainly there are other ways to devise 
integral preserving numerical schemes. A retraction is a map 

(j) : TM ^ M. 

Denote by (j)p the restriction of (j) to TpM and let Op be the zero- vector in TpM. Following [T], 
we impose the following conditions on (j) 

1. (j)p is smooth and defined in an open ball Br^iOp) c TpM of radius about Op. 

2. (I)p{v) = X if and only if = Op. 

3. ^Op^p = IdTpM- 

This implies in particular that ^p is a diffeomorphism from some neighbourhood ^ of Op to its 
image W = (j)p{^) c M. In what follows, we shall always assume that the step size used in 
the integration is sufficiently small such that both the initial and terminal point of the step are 
contained in such a set W. Furthermore, we shall assume that there is a given map c defined on 
some open subset of MxM containing all diagonal points (p^p)^ for which c{p^q) g M. Typically 
c(p, q) will be either p or or some kind of centre point between p and q to be defined later. We 
shall always require that c{p^p) = p for any p e M. 

We assume as before the existence of a bivector cu on M and a first integral H such that the 
differential equation can be written in the form ([2|. We introduce, for any pair of points p and 
q on M, an approximate bivector oj{p^q) such that 

Cd{p^p){v^w) = uj\p{v^w)^ \fv^weTpM. 

The discrete differential of a function H can now be defined for any pair of points (p^q) e M x M 
as a covector dH{p,q) e Tf^^ ^^M satisfying the relations 

H{q)-H{p) = {m{p,q),cl>-\<l)-re\p)) 
dH{p^p) = di^lp, for every p e M. 

where c = c{p^q) is the map referred to above. We define the integrator as 

x^+i = (/),(V^(x^,x^^i)), = ^-^x"") + /idi7(x^,x^^^) JcD(x^,x^^^) (16) 

It follows immediately that the method is symmetric if the following three conditions are satisfied: 

1. The map c is symmetric, i.e. c{p^q) = c{q,p) for all p and q. 

2. The discrete differential is symmetric in the sense that dH{p^q) = dH{q,p). 

3. The bivector uj is symmetric in p and q: oj{p^q) = uj{q^p). 
The condition [T]) can be achieved by solving the equation 

re\p) + 4>:\q) = 0, (17) 

with respect to c. 
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We can now write down a version of the AVF type discrete differential. Let 7^ = (1 - + 
where p = (l)c{v), q = (l)c{w). Then 

dff(p,g) = ^V:diJU(-,,)dC (18) 

Similarly, assuming that M is Riemannian, we can define the following counterpart to the Gon- 
zalez midpoint discrete gradient 

dff(p,,) = dF|. + ^Mz^Mz^^M,^ , = 4>-\q)-rJ{p)^nM. (19) 



where we may require that c{p^q) satisfies (17) for the method to be symmetric. 

Example 1. We consider the sphere M = S^~^ where we represent its points as vectors in 
of unit length, \\p\\2 = 1. The tangent space at p is then identified with the set of vectors in 
orthogonal to p with respect to the Euclidean inner product (•,•)• ^ retraction is 

Mv,)^^^ (20) 

its inverse is defined in the cone {q : (p, q) > 0} where 

Q _ 



A symmetric map c{p,q) satisfying (17) is simply 

c{p,q)-/^-. (21) 

the geodesic midpoint between p and q in terms of the standard Riemannian metric on S^~^ . We 
compute the tangent map of the retraction to be 



1 / (c + ix) (c + ii)\ 

l + u\\2 \ \\C + U\\l I 



As a toy problem, let us consider a mechanical system on S'^ . Since the angular momentum in 
body coordinates for the free rigid body is of constant length, we may assume (p^p) = 1 for all p 
and we can model the problem as a dynamical system on the sphere. But in addition to this, the 
energy of the body i preserved, i. e. 



which we may take as the first integral to be preserved. Here the inertia tensor isl= diagili ,12,13)- 
The system of differential equations can be written as follows 

^^^{dHAuj)\^^p^l-'p 

where the righthand side in both equations refer to the representation m R^. A symmetric con- 
sistent approximation to 00 would be 

^(P,^)(<^,/^) = X /^) 
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Figure 1: The free rigid body equation and its angular momentum in body coordinates. The 
trajectories are curves of constant energy on the sphere, computed by a method presented here. 
Moments of inertia used are I = diag(l, 2,3). 



We write = c + with the notation in (IS), this is a linear function of the scalar argument 
^. and thus, (pcile) ^ -^^/IK^II f'^om (20). We therefore derive for the AVF discrete gradient 

dH{p,q) = £ ^ (rVc(7«) - (^c(7«),rVc(7«))0c(7j)) 

This integral is somewhat complicated to solve analytically. Instead, we may consider the dis- 
crete gradient (19) where we take as Riemannian metric the standard Euclidean inner product 
restricted to the tangent bundle of S'^ . We obtain the following version of the discrete differential 
in the chosen representation 

dH{p,q) 



1 



m ■ 



\\Q-P\\ 



{H{q)-H{p)){q-p)], m = 



)• 



p^q 



The corresponding method is symmetric, thus of second order, and in Figure [7] we used this 
method to draw the trajectories of the free rigid body problem. 



4 Methods of higher order 

One viable way to obtain higher order variants of the proposed methods is by following the 
collocation strategy proposed by Hairer in [9 and Cohen and Hairer in [6], see also [2], and 
generalise it to the Lie group setting. 

Let ci, . . . , C5 be distinct real numbers such that < q < 1 and 6^ ^ for all i. We consider 
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the polynomial (j(^h) of degree s such that 

a(0) = 0, (22) 
dexp""^ {dHj -ic^j) , (Tj ■= a{cjh), (23) 



where 



dHj ■= (dexp^^^)' dexp;^^^) (Kxp(a(^/^))xo ^^exp(a(^/^))xo) (24) 

and cDj : g'' x ^ R are exterior 2-forms satisfying 

ujj{a,b) uJxARx-^ a,R\-ih), a,h e q\ (25) 

3 3 

with := exp(crj). The numerical solution after one step is xi = Rxq exp(cr(/i)). 
The collocation polynomial is obtained by integrating 

^Mih) = t ^jiOdeWa] {dHj ^Oj), (26) 

where ij j = 1 , . . . , s are the Lagrange basis functions and so 

.7 = 1 ^0 



Energy preservation 

To prove energy preservation of the proposed method we consider the path X{£^h) = Rxq exp{(j{£^h)) e 
G such that X(0) = xq and X{h) = Xi and integrate along this path obtaining 

using that 

— exp(a(^/i))xo = K^pi^a{^h))xo dexp^a/i) (^^(^^)) ' 



and (26), we obtain 



Jf-i ^ ^ 1 - 

{^^^Pa{^h)Rlxp(a(^h))xo ^^exp(cr(^/i))xo 5 




and further 

3 = 1 

SO that finally 



H{xi)-H{xo)=Y.^ji L -1^(^^^P^,^) ^e^Pa(^/i)Kxp(a(^/i))xo^^exp(a(^/i))xo^^.^^j ^^j). 
7 = 1 



H{xi)-H{xo) = Y,bj{dHj,dHjJCdj) = 0. 
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Order 



Considering the change of variables x{^h) = ex.p{Q{^h))xo ^ g [O^l]? and x the solution of 
([2|, by differentiation we get the following differential equation for Q: 



Depending on the choice of quadrature points and weights, the collocation method (22), (23) 
and ( 24 ) approximates ^7 at a certain order p, this suffices to guarantee that the overall method 
attains the same order, see also [6] and [M]. 

Example 2. Consider the collocation points ci^2 = | ^ ^ (the nodes of the Gauss quadrature) 
and let £i{(_) = {i - C2) I {ci - C2) and £2(0 = (^-Ci)/(c2 -Ci) be the corresponding Lagrange basis 
functions, let 

= fJh(r)dTdexp-l (dH^JOr) + £^ i2(r) dr dexp'] {^2^002) , 



and with dHj, j = 1,2^ given by (24), with hi = b2 = I, and Oj, j = 1,2^ given by (25). Then the 
one step method 

xi = exp(cr(/i)) Xq 
is a symmetric energy-preserving method of order 4. 



5 Examples and numerical experiments 

In the numerical experiments we consider the equations of the attitude rotation of a free rigid 
body in unit quaternions and a problem of elasticity the equations for pseudo-rigid bodies on 
the cotangent bundle of GL{3) (or SL{3)). 

5.1 Attitude of a free rigid body 

The set 

^^^{qGR^I ||q||2 = l}, q=[go,qf ,qeR^ 
with the quaternion product 

p • q [poqo - p^q, pop + ^oq + P x q]- 
is a Lie group. The corresponding Lie algebra is 

5^ :={[0,v]GR4|veR^}, 

and can be identified with R^. 

We consider the representation of the Euler attitude equations for the free rigid body in unit 
quaternions: 

q = /(q)-q, fin) = q-v-qc, 

and ^ 

V = [0,v], V = -r^f (qc)mo, 
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and I is the inertia tensor (a 3 x 3 fixed, diagonal matrix). The Euler- Rodriguez map S : ^ 
SO (3) is defined by 

f(q) :=/3 + 2^oq + 2q', 

where Is is the 3x3 identity matrix and q is defined as follows by means of the components of 

-qs q2 
Qs -qi 
-q2 qi 

The energy function preserved along q(t) is 

H(^) = imo£(q)r^f(qc)mo. 

We can choose the Euclidean inner product on := {[0, v] e | v g R^}, to play the role of the 
Riemannian metric on 5^, and the corresponding Riemannian gradient is 



gradi^ = {h - qq^) Vi^, 



with /4 the 4x4 identity. It follows that the bivector uo = Ig^f/^p can be expressed by the 4x4 
rank- 2 matrix 



^i?(q) 



II7IP 



where ^ = /(q) and 7 = gradJYl^ • qc. As a discrete bivector to construct the method we choose 
cj(q, q') = ujR{q), q = exp(7^/2)q, r] = log(q • q^). 
Identifying 5^ with its dual, and using the Gonzales trivialised discrete differential we get 

gradi^(q,q ) =7+ -— rj. 

mr 

The energy-preserving symmetric Lie group method is 



q' = exp (/icc;(q, q')gradi^(q, q')) • q. 

Alternatively we can use the averaged trivialised discrete differential obtained averaging the 
Riemannian gradient as 



gradi7(q,q0 = f^iOd^. liO = gradi^l^^^^ • q^CO, 



(27) 



and q(0 = exp(^ log(q' • qc)) • q. 

In figure [2] we apply a symmetric energy-preserving method of order 2 and 4 to the free rigid 
body problem in the formulation presented in this section, and we compare them with an explicit 
Lie group method of order 2 based on Heun's Runge-Kutta formula. We have used the discrete 
trivialized differential (27) in the Lie group method of order two and the formulation outlined in 
section [4] for the 4-th order variant of the method. The collocation points for this method are 
ci^2 ^ ^ ^ The integrals are approximated with accurate quadrature formulae. 
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Figure 2: Rigid body problem, qo = [1,0,0,0]^. Inertia tensor I = diag(l, 5, 60), mo = Ivq 
and vo = [1,0.5,-1]^. (Left) Order of the energy preserving methods, step-size versus norm of 
the global error: dashed-dotted line second order symmetric energy-preserving method based on 
the TDD ( [27| ); dotted line symmetric energy-preserving method of order 4 (section [4|; the solid 
lines are reference lines of order 2 and 4. (Right) momentum vector for the energy-preserving 
symmetric Lie group method of order 4, time interval [0,50], step-size h = . 



5.2 Pseudo-rigid bodies 

We consider the Hamiltonian equations describing rotating homogeneous elastic rigid bodies 
[12], [11]. A pseudo-rigid body is a three dimensional elastic body whose deformation gradient 
F = F{t) e GL^{3) is assumed to be constant throughout the body B: for all X e S c R^, the 
deformation is always given by the matrix vector product FX and F is not depending on X. 
In the incompressible case F e SL{3). The configuration space of a pseudo-rigid body is the 
group G equal to GL^{3) or S'L(3), and the fiow of the corresponding Hamiltonian equations 
evolves on T'^G ^ G x . With the semidirect-product group structure induced by the group 
multiplication in G, G x q"" is also a Lie group. In our particular example GL^{3) x Qi{3y , we 
identify Qi{3y with its dual, and we use coordinates (F^P) where F e GL^{3) and P € Qi{3) 
respectively. The corresponding Lie algebra is 0[(3) x ^[(3)^^. We denote with W(C) the stored 
energy function depending on the Cauchy-Green tensor C = F^ F, and with E the inertia tensor. 
The Hamiltonian function is 

H(F, P) := K{P) + H^(F^F), K{P) = ^ {P, PE-^), 

where (P^V) = tr(P-^V) is the standard matrix inner product, i.e. the duality pairing between 
between TG and T'^G. The canonical Hamilton's equations take the form 

where (ff , ff ) = (-2F VW(F^F), P£;-i). 

In our experiments we consider a St Venant-Kirchhoff material leading to a stored energy 
function 

W{C) = \ A(tr(C - I)f + Mtr((C - if), 

with tr denoting the trace operator and A and \i the Lame constants (such that /i > and 
3A + 2/i > 0, /i = 1 and A = | in the experiments). 
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The bivector uj is represented in this case simply by the 6x6 inverse Darboux matrix 



O I 
-I o 



The semi-direct product group multiphcation in G x g* is 

{Fi,Pi) ■ (F2, P2) = (F1F2, Pi + Ad^-i P2), 

and as matrix operation Ad^-i P2 = Pf ■^P2Pf-^, we denote with Exp and Log the exponential 
and logarithm between G x g* and g x g* respectively, these are defined by: 

Exp((^,/i)) = (exp(^),dexpj(/i)), Log({g,<j)) = (log(g), (dexp!iog(g) 

where exp and log are the corresponding maps for G. We consider 

V = iVF.Vp) = Log((Fi,Pi) • (Fo,Po)"'), 

and denote with (71,72) — ^(f P) ^^\(f p) ^ sK^)'^ ^ qK^) the trivialized differential in the 
point (F,P) = Exp(7/f/2,7/p)(Fo,^o), we have 



71 



SF 



-adl(P), 



72 



SP 



F- 



{F,P) 



and in matrix form ad^^ (^) - I2 ^ ~ ^72 • The trivialised discrete differential ( 15 ) becomes in 
coordinates 



\{F,P),{F',P')) = (71,72) +Q^(^F,^p), : 



H{F\P') - H{F,P) - (7i,7?p) - (7?p,72) 



and where the metric is deduced by the standard matrix inner product, 

\\r]f := tT{r]Fr]F)+tT{r]pr]p). 
The energy preserving method can then be formulated as 

(Pi, Pi) = Exp(/i(72 +a7?p,-7i -a7^F))(Po,Po). 



(30) 



In figure |3] we report the results of a simulation for this problem. We compare an explicit 
Lie group method (Heun's method), a symmetric Lie group method and the symmetric energy- 
preserving Lie group method presented in this section. The symmetric Lie group method is 



obtained by setting a = in (30). All methods are Lie group methods of order 2. We have halved 



the step-size for the explicit second order Lie group method, to avoid instability. All Lie group 
methods have the property that det(Pn) > for all n. The explicit Lie group method fails to 
preserve the energy of the problem (figure [3] top left), the symmetric Lie group methods have 
both a much smaller energy error (figure [3] top right), the symmetric energy preserving Lie group 
method preserves the energy to very high precision. We report here experiments with diagonal 
initial values for F and P. We have performed experiments also with non diagonal initial values 
obtaining similar results. The performance of the methods is relying on the accurate computation 
of matrix functions, and in particular matrix logarithms. In figure |3] (bottom right) we show the 
difference in the determinants of F for the two symmetric methods: the symmetric one denoted 
(sym) and the symmetric and energy preserving one denoted (EP). We plot det(P^^^)-det(P^^) 
for n = 1, . . . ,8000, this measures how the two numerical solutions depart from each other with 
time. 
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Figure 3: Pseudo-rigid body problem, Fq = Pq = diag(0. 2575, 0.8407, 0.2543), integration 
interval [0, 500]. Lame parameters A=|,/i = l. E = diag(l, 2, 3). (Top left) energy error explicit 
Lie group method (Heuns method), dotted line, h = 1/32, versus symmetric order 2 Lie group 
method, solid line, h = 1/16. (Top right) energy error symmetric order 2 Lie group method 
(energy error 10~^), versus the energy preserving Lie group method (energy error 10~^^, straight 
line), h = 1/16 for both methods. (Bottom left) deformation of the vector [1,1,1]^ under the 
transformation F(t), first 1000 steps, h = 1/16, energy-preserving method (with a different Pq). 
(Bottom right) difference of the determinants of F for the symmetric Lie group method and the 
symmetric, energy-preserving Lie group method: det(F^^^) - det(F^^), n = 1, . . . , 8000. 
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